#ifndef __ewpt__
#define __ewpt__
#include <iostream>
#include "TF1.h"
#include "TGraph.h"
#include "TMath.h"

using namespace std;


// a virtual class for storing information relating to potentialy
class EWPT{

 public:
  //declare constants
  const static double top_mass=173.3;
  const static double W_mass=80.39;
  const static double Z_mass=91.19;
  const static double h_mass=125.;
  const static double v=246.;

  //defined as mass*sqrt(2)/v
  static double top_couple;
  static double W_couple;
  static double Z_couple;
  static double h_couple;
  
  const static double PI=3.14159265358;
  //higgs quartic coupling

  //constants for high temperature expansion
  //bosonic:
  const static double b0 = -0.109662271123;
  const static double b2 = 1.0/24.0;
  const static double b3 = -0.026525823849;
  const static double b4 = 0.00856103393367;
  const static double b4log = -0.00316628698882;
  const static double b6 = 0.00001606809223562;
  
  //fermionic:
  const static double f0 = -0.09595448723281;
  const static double f2 = 1.0/48.0;
  const static double f4 = -0.00417162813538;
  const static double f4log = 0.0031662869888230;
  const static double f6 = -0.00011247664564935;

  EWPT();

  double Tc() const;
  //compute the critical temperature

  

  TGraph* GetGraph(double T=0) const;
  //get a potential given a temperature

  TGraph* GetIntGraph(int spin=0) const;
  //Get a graph of the integral, for debugging

  TGraph* GetLowGraph(int spin=0) const;
  //Get a graph of the integral, for debugging

  TGraph* GetHighGraph(int spin=0) const;
  //Get a graph of the integral, for debugging
  
  // protected:
  
  TF1* V_tree;
  TF1* V_loop_top;
  TF1* V_loop_W;
  TF1* V_loop_Z;
  TF1* V_loop_H;

  const static double matching_high=1.2;
  const static double matching_low=1.5;  

  
  virtual double operator()(double x, double T=0) const =0;
  //return zero temperature coleman weinberg potential
  
  double SM(double x, double T=0) const;
  //return the SM potential

  double Dx(double x, double T, double step=1e-3) const;
  //dx (partial derivative) of the potential
  //computed using discretization

  double Dx2(double x, double T, double step=1e-3) const;
  //dx2 (second partial derivative) of the potential
  //computed using discretization

  double find_v(double T, double step=1e-3, double guess=500.)
    const;
  //find the nonzero local minimum

  double Tc(double& vev) const;
  //compute the critical temperature and return the critical vev

  void Tc_helper(double& T, double& vev, double step) const;
  //compute the critical temperature with some acuracy
  //helper function

  double tint_boson(double r) const;
  //return the thermal integral for bosons
  //i.e. int dx x^2 log(1-e^{x^2+y^2})
  //input = y = mh^2/T^2

  double tint_fermion(double r) const;
  //return the thermal integral for fermions
  //i.e. int dx x^2 log(1-e^{x^2+y^2})
  //input = y = mh^2/T^2

  //sums of bessel functions for approximation
  double bessel_sum(double r) const
  {
    return r*r*TMath::BesselK0(r) + 2*r*TMath::BesselK1(r);
  }

  double smooth(double r) const;
  //return a double between 0, 1. 
  //Used for smoothingly interpolating two func


};

#endif
